Canal	Longitud de onda Central \(\mu m\)	Observaciones	Función de peso
7	3.9	Sensible tanto a la radiación del sol como a la emitida por la Tierra	Maxímo en superficie
8	6.2	Sensible al vapor de agua	Máximo en 425 hPa.
9	6.9	Predomina la absorción por vapor de agua	Maxímo cercano a 460 hPa.
10	7.3	Sensible al vapor de agua	Máximo en 640 hPa.
11	8.4	Para estudio de microfísica de nubes y plumas volcánicas con ceniza y dióxido de sulfuro	Máximo en niveles bajos
12	9.6	Absorción por ozono	Máximo en la estratósfera
13	10.3	Ventana radiativa, absorción despresiable en la atmósfera. Cerca del \(\lambda\) de máxima emisión terrestre	Máximo en superficie
14	11.2	Ventana radiativa con mayor absorción de vapor de agua
15	12.3	Ventana radiativa con mayor absorción de vapor de agua
16	13.3	Absorción por dióxido de carbono	Máximo cerca de superficie

Corrección del BIAS

De acuerdo a la guía del usuario avanzda de GSI, para asimilar radiansas y que estas generen un impacto positivo en el prónosticos, es impresindible corregigir el bias asociado a cada instrumento y canal.

GSI corrige el BIAS sobre la diferencia entre la observación y el background teniendo en cuenta dos posibles fuentes de BIAS:

Air mass bias correction (corrección variacional): usa los coeficientes definidos en satbias_in
Angle dependendat bias correction: asociado al ángulo de escaneo para cada punto y cuyos coeficientes se guardan en satbias_angle

Ambos set de coeficientes se deben ir actualizando en cada ciclo de asimilación, actualmente esa actualización se hace internamente (siempre que se configuren los argumentos necesarios en el namelist). De acuerdo a la guía los coeficientes para la correcciónn del bias tienen un spin up muy lento y deben ser entranados por semanas o meses. Una alternativa posibles es repetir cada ciclo de asimilación 10-12 veces antes de pasar al siguiente, cada vez actualizando los coeficientes.

Una tercera opción sería probar con los coeficientes que se guardas en GDAS y que usa GFS. No necesariamente es una buena solución ya que son modelos distintos con distinta resolución y dominio. Pero tal vez se puede partir de esos coeficientes para que el spin up sea más corto (en comparación con empezar en 0).

Queda por responder ¿cómo se que los coeficientes ya están bien entrenados y puedo usarlos?

Funciones de peso: http://tropic.ssec.wisc.edu/real-time/amsu/explanation.html

Lectura de coeficientes de GDAS

Maravillosamente los archivos en GDAS son de texto plano y tienen el formanto GSI friendly. Ahora me resulta obvio, pero nunca se sabe!

Entonces, hay un archivo cada 6 horas que se va actualizando.

gdas.txxz.abias y gdas.txxz.abias_pc: Según la guía de GSI es “combined satellite angle dependent and mass bias correction coefficient file”. El segundo para canales “pasivos” (que se monitorean).

Al parecer y de acuerdo a las versiones recientes de la guia de usuario, la correción del bias asociado al ángulo y el ¿modelos? ¿backgroud? se hace conjuntamente y por eso los coeficientes se combinan en un solo archivos.

# Es un archivo de texto plano, que emoción!!
satbias_in_file <- "/home/paola.corrales/datosmunin2/nomads.ncdc.noaa.gov/GDAS/201811/20181120/gdas.t00z.abias"

Para que la corrección se haga en conjunto es necesario indicarlo en el namelist: adp_anglebc = .true.

Archivos presentes en la salida de GSI

satbias_in: the input coefficients of bias correction for satellite radiance observations.
satbias_out: the output coefficients of bias correction for satellite radiance observations after the GSI run.
satbias_pc: the input coefficients of bias correction for passive satellite radiance observations
satbias_pc.out

Prueba de asimilación

ana_file <- "/home/paola.corrales/datosmunin/EXP/analysis.ensmean"
ana <- ReadNetCDF(ana_file, vars = c(p = "P", "PB", t = "T", qv = "QVAPOR", 
                                     lon = "XLONG", lat = "XLAT")) %>%
  .[, p := p + PB] %>%
  .[, t := tk(t, p)] %>%
  .[, rh := rh(qv, p, t)] %>% 
  .[, td := td(qv, p) + 273.15] %>% 
  .[, ":="(Time = NULL,
           west_east = NULL,
           south_north = NULL,
           qv = NULL,
           PB = NULL)] %>% 
  .[, c("u", "v") := uvmet(ana_file)] %>% 
  .[, file := "ana"]

guess_file <- "/home/paola.corrales/datosmunin/EXP/wrfarw.ensmean"
guess <- ReadNetCDF(guess_file, vars = c(p = "P", "PB", t = "T", qv = "QVAPOR", 
                                         lon = "XLONG", lat = "XLAT")) %>%
  .[, p := p + PB] %>%
  .[, t := tk(t, p)] %>%
  .[, rh := rh(qv, p, t)] %>% 
  .[, td := td(qv, p) + 273.15] %>% 
  .[, ":="(Time = NULL,
           west_east = NULL,
           south_north = NULL,
           qv = NULL,
           PB = NULL)] %>% 
  .[, c("u", "v") := uvmet(ana_file)] %>% 
  .[, file := "guess"]

dt <- rbind(ana, guess)

dt[bottom_top %in% seq(1, 35, 4)] %>% 
  dcast(bottom_top + lon + lat ~ file, value.var = "t") %>% 
  .[, diff_t := ana - guess] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = diff_t)) +
  scale_color_divergent() +
  geom_mapa() +
  coord_sf(xlim = c(-78, -50), ylim = c(-42, -19)) +
  facet_wrap(~ bottom_top) +
  theme_minimal()

dt[bottom_top %in% seq(1, 35, 4)] %>% 
  dcast(bottom_top + lon + lat ~ file, value.var = "p") %>% 
  .[, diff_p := ana - guess] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = diff_p)) +
  scale_color_divergent() +
  geom_mapa() +
  coord_sf(xlim = c(-78, -50), ylim = c(-42, -19)) +
  facet_wrap(~ bottom_top) +
  theme_minimal()

Actualización de coeficientes

En la etapa del LETKF la subrutina update_biascorr() se encarga de actualizar los coeficientes que se usan en la corrección del bias aplicando la metodología propuesta por Miyoshi et. al. (2010) y que según Zhu et. al. (2019) en análoga a la versión variacional que utiliza GSI. En el mismo paper indican que la optimización con filtro de kalman en algunos casos funciona mejor y esto podría doberse a que es menos sensible al tamaño de la muestra (cantidad de observaciones).

La rutina calcula el incremento de los coeficientes \(\delta \beta\) de la siguiente manera:

\[\delta \beta = (B_{\beta}^{-1} + P R^{-1}P^{T})^{-1}PR^{-1}[y-H(x)-P^T\beta]\]

Donde \(B^{-1}_\beta\) la matriz de covarianza de los coeficientes es diagonal y sus elementos valen 0.1.

La estimación de \(\delta \beta\) se realiza iterativamente tantas veces como el argumento numiter del namelist lo indique (por defecto 6).

Previamente al cálculo del incremento se realiza una estimación adaptativa de la varianza de los errores que impacta en la matriz \(B\). Actualmente el algoritmo tiene en cuenta:

Si es la primera vez que se asimila ese canal/sensor (revisar): \(biaservar = 10000\)
Si en algún ciclo no hay observaciones de un canal/sensor: \(biaservar = 2*biaservar + 1e-6\) (con un tope de \(biaservar = 10\))
O, si eso está desectivado puede tomar un valor fijo dividido la cantidad de observaciones.

Predictores

Según el código los predictores \(p_i\) son:

pred(1,:) = global offset
pred(2,:) = zenith angle predictor, is not used and set to zero now
pred(3,:) = cloud liquid water predictor for clear-sky microwave radiance assimilation
pred(4,:) = square of temperature laps rate predictor
pred(5,:) = temperature laps rate predictor
pred(6,:) = cosinusoidal predictor for SSMI/S ascending/descending bias
pred(7,:) = sinusoidal predictor for SSMI/S
pred(8,:) = emissivity sensitivity predictor for land/sea differences
pred(9,:) = fourth order polynomial of angle bias correction
pred(10,:) = third order polynomial of angle bias correction
pred(11,:) = second order polynomial of angle bias correction
pred(12,:) = first order polynomial of angle bias correction

Se calculan utilizando información del guess y de la observación de radianza y luego se utilizan para calcular la corrección del bias para cada observación y que será aplicada sobre el incremento.

\[ tbc = tbcnob + \sum_{i} p_i (\beta + \delta \beta)\]

files <- list.files("analisis/diagfiles/E7", pattern = "asim", full.names = TRUE)
diag <- map(files, function(f){
  meta <- unglue::unglue(f, "analisis/diagfiles/{exp}/asim_{sensor}_{plat}_{date}.ensmean")
  #print(f)
  out <- fread(f)
  # .[V10 == 1] %>% 
  
  if (file.size(f) != 0) {
    out[, date := ymd_hms(meta[[1]][["date"]])]
  }
  out
}) %>%
  rbindlist()

## Warning in fread(f): Stopped early on line 276. Expected 28 fields but found
## 27. Consider fill=TRUE and comment.char=. First discarded non-empty line:
## <<mhs_n19 1 -23.43 304.42 455.39 -0.09 99999997952.000000 99999997952.000000
## 99999997952.000000 0.000000 50.000000 0.950000 1.974102 -46.150002 82.910004
## 1.000000 0.000000 0.000000 0.000000******************** 0.545147E+00 0.100000E
## +01 0.000000E+00 0.000000E+00 0.169061E+01 0.130024E+01 0.000000E+00 0.000000E
## +00>>

## Warning in fread(f): Stopped early on line 896. Expected 28 fields but found
## 27. Consider fill=TRUE and comment.char=. First discarded non-empty line:
## <<mhs_n19 1 -22.45 287.04 0.00 -0.40 99999997952.000000 99999997952.000000
## 99999997952.000000 0.000000 50.000000 0.591223 0.283433 6.940000 259.660004
## 0.000000 0.000000 0.000000 9.000000******************** -0.155139E+00 0.100000E
## +01 0.000000E+00 0.000000E+00 0.152438E+00 -0.390433E+00 0.000000E+00 0.000000E
## +00>>

colnames(diag) <- c("sensor", "channel", "lat", "lon", "press", "dhr", "tb_obs", "tbc", "tbcnob",
                    "errinv", "qc", "emis", "tlapchn", "rzen", "razi", "rlnd", "rice", "rsnw", "rcld", 
                    "rcldp", paste0("pred", seq(8)), "date")

diag[sensor == "amsua_n15" & channel %in% c(1:10)] %>% 
  melt(measure.vars = c("tbcnob", "tbc")) %>% 
  .[, .(mean = mean(value),
           sd = sd(value)), by = .(variable, channel, date)] %>% 
  melt(measure.vars = c("mean", "sd"), variable.name = "estadistico") %>% 
  ggplot(aes(date, value)) +
  geom_point(aes(color = variable)) +
  geom_path(aes(color = variable, linetype = estadistico)) +
  scale_color_brewer(name = NULL, palette = "Set1", labels = c("tbcnob" = "sin BC",
                                                  "tbc" = "con BC")) +
  # geom_vline(xintercept = 0) +
  facet_wrap(~channel, scales = "free") +
  labs(x = "Obs - Guess") +
  theme_minimal() +
  theme(legend.position = "bottom")

diag[sensor == "amsua_n15" & channel %in% c(1:8)] %>% 
  .[date == date[1]] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_point(aes(color = tbc - tbcnob)) +
  scale_color_divergent(name = "BIAS") +
  geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  labs(title = "amsua_n15",
       x = "", y = "") +
  facet_wrap(~channel, ncol = 4) +
  theme_minimal()

Pred2 = 1 Pred3 = 0 Pred7 = 0 Pred8 = 0

diag[sensor %in% c("amsua_n15")] %>% 
  ggplot(aes(factor(channel), tbc - tbcnob)) +
  geom_hline(yintercept = 0, color = "darkorange") +
  geom_boxplot() +
  labs(x = "Channel", y = "BIAS") +
  facet_wrap(~sensor) +
  theme_minimal()

# diag[sensor %in% c("airs_aqua") & channel < 300] %>% 
#   ggplot(aes(channel, tbc - tbcnob)) +
#   geom_boxplot(aes(group = channel))

Prueba de diagfiles en formato ncdf

files <- list.files("analisis/diagfiles/", pattern = "amsua_n15", full.names = TRUE)

diag <- map(files, function(f) {
  var_names <- GlanceNetCDF(f) %>% 
    .$vars %>% 
    names()
  
  diag <- ReadNetCDF(f, vars = var_names[c(10, 12, 13, 55:79)]) %>% 
    .[, file := f]
}) %>% 
  rbindlist()

diag[Channel_Index %in% c(1:4, 6, 7)] %>% 
  melt(measure.vars = c("Obs_Minus_Forecast_unadjusted", "Obs_Minus_Forecast_adjusted")) %>% 
  ggplot(aes(value)) +
  geom_density(aes(color = variable)) +
  geom_vline(xintercept = 0) +
  facet_wrap(~Channel_Index, scales = "free") +
  labs(x = "Obs - Guess") +
  theme_minimal() +
  theme(legend.position = "bottom")

diag %>%
  melt(measure.vars = patterns("^BCPred_")) %>%
  ggplot(aes(ConvertLongitude(Longitude), Latitude)) +
  geom_point(aes(color = value)) +
  scale_color_viridis_c() +
  geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  facet_wrap(~variable, ncol = 4) +
  theme_minimal()

diag[Channel_Index == 2] %>% 
  
  ggplot(aes(ConvertLongitude(Longitude), Latitude)) +
  geom_point(aes(color = Observation)) +
  scale_color_viridis_c() +
  geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  # facet_wrap(~variable, ncol = 4) +
  theme_minimal()

104_radiancias